pygsf 1: spatial data

March-April, 2018, Mauro Alberti, alberti.m65@gmail.com

Developement code:



In [3]:

    
%load_ext autoreload
%autoreload 1









    



The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

1. Introduction

gsf is a library for the processing of geometric and geographic data, with a focus on structural geology data. It is composed by a Python 3 module, pygsf, and an experimental, still in progress Haskell module, named hsgsf. This notebook will present the Python version.

Since we will plot geometric data into stereonets, prior to any other operation, we import mplstereonet and run the IPython command %matplotlib inline, that allows to incorporate Matplotlib plots in the Jupyter notebook.



In [4]:

    
%matplotlib inline

We import all classes/methods from the geometry sub-module:



In [5]:

    
from pygsf.spatial.vectorial.vectorial import *

2. Basic spatial data types: points and planes

The reference axis orientations used in pygsf are the x axis parallel to East, y parallel to North and z vertical, upward-directed.

2.1 Cartesian points

A point is created by providing three Cartesian coordinates:



In [6]:

    
p1 = Point(1.0, 2.4, 0.2)  # definition of a Point instance

We calculate its distance from the reference frame origin:



In [7]:

    
p1.len3D  # distance of a point from the origin









    Out[7]:





2.6076809620810595

When considering two points, we can calculate their 3D distance as well as their horizontal, 2D distance:



In [8]:

    
p2 = Point(0.9, 4.2, 10.5)



In [9]:

    
p1.dist3DWith(p2)  # 3D distance between two points









    Out[9]:





10.45657687773585



In [10]:

    
p1.dist2DWith(p2)  # horizontal (2D) distance between two points









    Out[10]:





1.8027756377319948

Among other methods, we can:

translate the point position by providing three offset cartesian values (x, y and z) or directly via a vector;
check if two points are within a given range of each other;
convert a point to a vector.

2.2 Cartesian planes

A Cartesian plane can be defined in a few different ways, with the simplest one by providing three points within the plane:



In [11]:

    
pl1 = CPlane.fromPoints(Point(0, 0, 0), Point(1, 0, 0), Point(0, 1, 0))



In [12]:

    
print(pl1)









    



CPlane(0.0000, 0.0000, 1.0000, 0.0000)

The four coefficient returned (a, b, c and d) define the Cartesian plane as in the equation:

ax + by + cz = d

For the given example, the equation is satisfied for all x and y values provided z is zero. We are therefore considering a horizontal plane passing through the frame origin.

The versor normal to a Cartesian plane is obtained by the method:



In [13]:

    
normal_versor = pl1.normVersor()  # versor (unit vector) normal to the provided Cartesian plane



In [14]:

    
print(normal_versor)









    



Vect(0.0000, 0.0000, 1.0000)

In this example the normal versor is vertical.

We can calculate the intersection, expressed as a versor, between two Cartesian planes:



In [15]:

    
pl1, pl2 = CPlane(1, 0, 0, 0), CPlane(0, 0, 1, 0)
inters_v = pl1.intersVersor(pl2)  # intersection versor between two Cartesian planes 
print(inters_v)









    



Vect(0.0000, -1.0000, 0.0000)